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Abstract. We propose an improvement of the differential method for the computation of the equa- 
tion of state of QCD from lattice simulations. In contrast to the earlier differential method our 
technique yields positive pressure for all temperatures including in the transition region. Employing 
it on temporal lattices of 8, 10 and 12 sites and by extrapolating to zero lattice spacing we obtained 
the pressure, energy density, entropy density, specific heat and speed of sound in quenched QCD for 

0. 9 < T/T c < 3. A comparison of our results is made with those from the dimensional reduction 
approach and a conformal symmetric theory at high-temperature. 
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1. Introduction 

There is growing acceptance of the view that in the ongoing experiments in Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven a new form of matter has been created [1]. This 
new form of matter is thought to be a fluid of strongly interacting quarks and gluons. In 
lattice studies of quenched QCD it was found earlier that the entropy density s [2,3] and 
the mean free time r, derived from the electrical conductivity [4], together gave rise to 
a dimensionless number rs 1 / 3 0.8 [5]. In the non-relativistic limit this dimensionless 
number measures the mean free path in units of interparticle spacing, and is therefore large 
in a gas but of order unity in a liquid. This indicated that the deviation of the energy density 
(e) and pressure (P) in the high temperature phase of QCD from their ideal gas values may 
be due to a previously underappreciated feature of the plasma phase — that it is far from 
being a weakly interacting gas. 
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Earlier expectations that a weakly interacting gas of quarks and gluons would be formed 
in the experiments were based on perturbative calculations [6] which failed to reproduce 
these lattice results [2]. There have been many suggestions for the physics implied by the 
lattice data — the inclusion of various quasi-particles [7], the necessity of large resumma- 
tions [8], and effective models [9] being a few. Investigation of screening masses also gave 
evidence for strong departure from perturbative results [10-14]. Interestingly, there has 
been a suggestion that conformal field theory comes closer to the lattice result [17]. This 
assumes more significance in view of the fact that a bound on the ratio of the shear vis- 
cosity and the entropy density, s, conjectured from the AdS/CFT correspondence [18] lies 
close to that inferred from analysis of RHIC data [19] and its direct lattice measurement 
[20,21] as well as the lattice results of a different transport coefficient [4]. 

The equation of state (EOS) is one of the most basic inputs into the analysis of experi- 
mental data. Two decades ago, a method was devised to compute the EOS of QCD on the 
lattice [22]. However, soon it was found [23] that this method yielded negative P near the 
critical temperature, T c . At that time it was thought that this problem of the "differential 
method", as it is called now, is solely due to the use of perturbative formulae for various 
derivatives of the coupling. To cure this problem of negative pressure the non-perturbative 
"integral method" was introduced [24,2]. It bypasses the use of perturbative couplings 
by employing the thermodynamic relation F = —PV and using a non-perturbative but 
phenomenologically fitted QCD /3-function. If the EOS were to be evaluated by the inte- 
gral method then fluctuation measures (e.g. the specific heat at constant volume C v ) can 
only be evaluated through numerical differentiation, which is prone to large errors [25]. 
Moreover, the relation F — —PV assumes the system to be homogeneous. Since the pure 
gauge phase transition in QCD is of first order the system is not homogeneous at T c . Thus 
one makes an unknown systematic error in the integral method computation by integrating 
through T c . This is in addition to a small systematic error due to setting P = just below 
T c and the numerical integration errors. Clearly, our confidence in the lattice results on 
the EOS would be boosted if an entirely different method of EOS determination yields the 
same results: it would tantamount to a good control over many systematic errors in both. 

In this paper we propose a modification of the differential method which gives positive 
pressure over the entire temperature range for even relatively coarser lattices . We choose 
the temporal lattice spacing (o T ) to set the scale of the theory, in contrast to the choice of the 
spatial lattice spacing (a s ) in the approach of [22]. This change of scale is analogous to the 
use of different renormalization schemes. As a consequence, our method could be called 
the t-favoured scheme and the method of Ref. [22] may be called the s-favoured scheme. 
In fact, in a different context, this choice of scale has already been used in Ref. [26]. Here 
we show that this choice leads to positive pressure for the entire temperature range, even 
when one uses one-loop order perturbative couplings. Since the operator expressions are 
derived with an asymmetry between the two lattice spacings a s and a T , the s-favoured and 
t-favoured schemes give different expressions for the pressure. In that sense the use of 
t-favoured scheme is tantamount to the use of better operators. 

Being a differential method the t-favoured scheme can be easily extended for the cal- 
culation of fluctuation measures like C v , following the formalism developed in Ref. [3]. 
In a theory with only gluons there is only this one fluctuation measure. Related to this 
is a kinetic variable, the speed of sound, C s , which can also be evaluated in any operator 
method. We report measurements of both in the temperature range 0.9T C < T < 3T C 
through a continuum extrapolation of results obtained using successively finer lattices. 
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Not only do these quantities provide further tests of all the models which try to explain 
the lattice data on the EOS they also have direct physical relevance to experiments at RHIC. 
In a canonical ensemble the specific heat at constant volume is a measure of energy fluc- 
tuations. It was suggested in Ref. [27] that event-by-event temperature fluctuation in the 
heavy-ion collision experiments can be used to measure C v . The speed of sound, on the 
other hand, controls the expansion rate of the fire-ball produced in the heavy-ion collisions. 
Thus the value of C s is an important parameter in the hydrodynamic studies. It has been 
noted that the magnitude of elliptic flow in heavy-ion collisions is sensitive to the value of 
C s [28]. 

The measurement of C v and C s also directly test the relevance of conformal symmetry 
to finite temperature QCD. QCD is known to generate the scale, Aq C d, dynamically and 
thus break conformal invariance. The strength of the breaking of this symmetry at any scale 
is parametrized by the /3-function. An effective theory which reproduces the results of 
thermal QCD at long-distance scales could still be close to a conformal theory. The result 
of Ref. [17] for the entropy density, s, in a Yang-Mills theory with four supersymmetry 
charges (Af = 4 SYM) and large number of colours, N c , at strong coupling, is 

- = .f(g 2 N c ), where s = ^r 2 A 2 T 3 and 
s 3 

m - \ + |c(3)z- 3/2 + • • • , a) 

g being the Yang-Mills coupling 1 . For our case of N c — 3, the well-known result for the 
ideal gas, s = 4(A 2 - l)7r 2 T 3 /45 takes into account, through the factor - 1, the 
relatively important difference between a SU (N c ) and an U (N c ) theory. 

The paper is organized as follows. In the next section we present the formalism and lead 
up to the measurement of C v and C 2 on the lattice in Section 2.2. In Section 3 we give 
details of our simulations and our results. Finally, in Section 4 we present a discussion of 
the results. 



2. Formalism 

Various derivatives of the partition function, Z(V, T), where V is the volume and T the 
temperature, lead to thermodynamic quantities of interest. In particular the energy density, 
e, and the pressure, P, are given by the first derivatives of In Z, 



T\ T d\nZ(V,T) 



V dT 



and P=[ r\v dx * z{ y^ 



V dV 



(2) 



The second derivatives are measures of fluctuations. In the absence of chemical potentials 
a change of volume of a relativistic gas alters its pressure by changing particle numbers. As 
a result there is only one second derivative, namely the specific heat at constant volume — 



1- We thank Igor Klebanov for pointing out that the factor x 3/2 in the right hand side of Eq. (1) 
appears in early literature as (2x)~ 3 ^ 2 due to a different normalization of g 2 N c . 
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C v 



dT 



(3) 



Using thermodynamic identities, the expression for the speed of sound can be recast in 
the form 



dP 


dP 




f de 




Ik 


-dT 


J 


y dT 


J 



c 



where we have used the thermodynamic identity 



dP 
dT 



dS_ 

dV 



and dV 



s/T 3 

Cy/T^ 



e + P 



(4) 



(5) 



in conjunction with the definition of the total entropy, S, and the entropy density, s, above. 
Note that all these relations are valid for full QCD with dynamical quarks (without quark 
chemical potentials) as well as in the quenched approximation which this work deals with 
exclusively. 

A caveat about the first equality in Eq. (4) is in order. This remarkable formula (a 
generalization of a result first obtained in 1687 by Newton) equating a kinetic quantity, 
Cg, to a thermodynamic derivative is true for a homogeneous system. For a phase mixture 
at a first order phase transition there are kinetic processes, such as condensation of a fog, 
which cause this formula to break down [29] . The lore that C 2 S = at T c is due to the overly 
naive argument that P remains continuous while e undergoes a discontinuous change. In 
fact, the best that thermodynamics can do is to evaluate this formula in a limiting sense as 
one approaches T c either from above or below. The values of C s in these two limits need 
not even be continuous at a first order transition [30]. 



2.1 Energy density and pressure 



In order to distinguish between T and V derivatives, the differential method formulate 
the theory on a d + 1 dimensional asymmetric lattice having different lattice spacings in 
the spatial (a s ) and the temporal (a r ) directions. If the number of lattice sites in the two 
directions are N s and N T , then T = (A^T-a,-) -1 is the temperature and V = (N s a s ) d is the 
volume of the system. The derivatives needed for the thermodynamics are — 



d_ 

dT 



d 
da T 



and V 



_d_ 

dV 



Os_ _d_ 

d da s 



(6) 



holding N s and N T fixed. 

In the t-favoured scheme we introduce the anisotropy parameter £ and the scale a by the 
relations, 



— , and 



a = a T . 



(7) 



The partial derivatives with respect to T and V can then be written in terms of these new 
variables as 



4 



Pramana - J. Phys., Vol. xx, No. x, xxx xxxx 



Lattice QCD equation of state 



dT 



d_ 



d_ 

da 



and 



V 



dV 



d d£ 



(8) 



One obtains the second expression by writing a s = a£ and taking a partial derivative 
keeping a fixed. For the first expression, one takes a derivative with respect to a and then 
introduces constraints on the differentials d£ and da in order to keep a s fixed. This choice 
of scale a = a T seems to be natural, since most numerical work at finite temperature sets 
the scale by T = l/N T a T . For example, continuum limits are taken at fixed physics by 
keeping T fixed while changing N T and a T simultaneously. This is done not only when 
symmetric lattices are used, but also when the simulation is performed with asymmetric 
lattices [31]. 

In the s-favoured method [22], by contrast, the scale of the theory is set by the spatial 
lattice spacing, a = a s , at every £ and only after taking the £ — > 1 limit does the natural 
choice of scale emerge. The corresponding derivatives in this case are 



d 

T t— 

dT 



and 



d 




£ d 


a d 


dV 


T 


d <9£ 


a d da 



(9) 



On the anisotropic lattice the partition function of a pure gauge SU (N c ) theory with the 
Wilson action is defined as 



Z(V,T) = j VUe- s M, where 

d d 

S[U]=K S ]T Pij(x)+K T p oi(x)- 

x,ij—l x.i—1 



(10) 



Periodic boundary conditions are imposed in all directions. The plaquette variables are 
Pap(x) = 1 — Re trU a p(x), U a [}(x) is the ordered product of link matrices taken 
anticlockwise around the plaquette, starting at the site x and in the plane specified by 
the directions a and (3. We introduce the notation for the average plaquettes P s = 
2j2Pij(x)/d(d - 1)N*N T and P T = Y,Pnt{x)/dNfN T . Since the plaquette opera- 
tors have no explicit dependence on a and £ the derivatives with respect to these quantities 
vanish. The couplings may be written as 



leading to 



K, = 



2N C 



and 



-K. + 2N, 



9g; 



d£ 



2N C £ 



and 



(11) 



dK T 



K T + 2N C (- 



(12) 



Next, using the derivatives in Eq. (8) along with the definitions of P and e (see Eq. 2) 
one obtains, from the partition function of Eq. (10), the expressions 



a d+1 e 



d_ 
1 



i 



2 

d- 1 



£K' s D a + tK' T D T 
^K' S D S + iK' T D T 



+ 



d_ 



1 dK dK T 
■ a— — D s + a— — D T 



da 



da 



and 
(13) 
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where primes denote derivative with respect to £. In order to remove the trivial ul- 
traviolet divergence in these quantities, present even in the free case, a subtraction of 
the corresponding T = values is made, yielding Di — (Pi) — (Pq) above. Here 
P = 2^ P a p(x) I 'd(d + l)NfN T is the average plaquette value at T = 0, evaluated 
with periodic boundary conditions in all directions and with very large N T = N s . 
To determine the couplings K[ we use the weak coupling definitions [32] 

^ = m +Cii0+0[92{a)] ( ^ s ' r) - (14) 

With the condition that Cj(£ = 1) = 0, this is actually an expansion of the anisotropic lat- 
tice couplings gi(a, £) around the isotropic lattice coupling g(a). With the usual definition, 
a s = g 2 /47T, the /3-function is — 

d/ n M da s . . dg- 2 B(a s ) 

B{as)= 2^ glvmg a ^T = ^T (15) 

For a 3 + 1 dimensional theory one has B(a s ) = —(33 — 2Nf)a 2 /12ir + O(a^). In terms 
of the functions c s and c T introduced in Eq. (14) and the /3-function above one can rewrite 
the derivatives of the couplings as 



8K S N c B(a s ) 8K S , 

a ^r~ = b - ^' and = ~ K s + 2N c c s, 

da ^ajt, ot, 



8K T N£B{a s ) .dK T 2 , 

a—— = , and ^— — = K T + 2N c ^c T . (16) 

oa iraj at, 

The quantities c' s and c' T have been computed to one-loop order in the weak coupling limit 
for SU(N C ) gauge theories in 3+1 dimensions [33]. 



2.2 The specific heat and speed of sound 



It was pointed out in Ref. [3] that the specific heat can be most easily obtained by working 
with the conformal measure, 



C = 



and 



ac 

df 



(17) 



where A = e — 3P. Then, using Eqs. (4, 5, 17) it is straightforward to see that 



CV = / e/T d+l 
T d ~ \ p/T d+1 

f P/T d+1 
V e/T d +! 



C 



s T e 
+ 'dT d + 1 
Te / T d+i- 



and 



ds/T d 



(18) 



One needs the expression for V in terms of the plaquettes in order to proceed. To this 
end we introduce the two functions 
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F(Z,a) = ^-J- = a 



G(H,a) 



d 

-ea d+1 t d 
d 



d-ldK, 8K T 



2 da 
d-l 



da 

K'D S + K'D T 



and 
Ffoa). 



Since C = —F/G, one finds that 



L F 8T 



T8G 
+ L G 8T 



(19) 



(20) 



The derivatives of F and G will involve the variances and covariances of the plaquettes 
and the second derivatives of the couplings. These second derivatives of the couplings 
are — 



,d 2 K s 
da 2 



da 

a djK' T _ CB(a s ) 
da 2na 2 
B(a s ) dK s 



2na 2 ^ 



da ' 



, d 2 K T £B(q 8 ) = dK T 
da 2 2na 2 da 



The numerical values of c" 's have been evaluated in Ref. [3]. 

Turning now to the derivatives of F and G in Eq. (19) one obtains 



c dF 

dF 
da 



+ a' 



d-18K>. ^ 



2 da 
d-!8K s 
~ da~ 



da 



d-\dKs D , | dK TD , 



2 da 
d - 1 d 2 K t 
~~ 2 IkP 



da 



d - 1 dK s dD s dK T dD T 



2 da da 
Also from Eq. (19) it follows 
'd-l 



da da 



d£, 

+e 

dG 

a— = £a 
da 

+ & 



2 

d-l 



K' S D S + K' T D T 
K' S D' S + K' T D' T 



d-l 



JF 



k' s 'd s + k*;d t 



and 



d-ldK' dK' T 
da 

,dD T 



da 



da 



da 



dF 
da 



(21) 



and 



(22) 



(23) 



Since the plaquette operators do not explicitly depend on £ and a one can easily take the 
derivatives of the vacuum subtracted plaquette expectation values. These are 
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ZD'i = -dN T N d s 

3D; 



and 



da 



-dN T N" 



d-l dK s dK T 
a— — a si + a— — a T 



da 



da 



(24) 



where cry = (DiDj) — (Di)(Dj). Throughout this paper we will refer to try (i ^ j) as 
'variances of plaquettes' and an as 'covariances of plaquettes'. Note that Eq. (18) implies 
that C v and C s should be independent of the volume. Consistent with this, the derivatives 
in Eqs. (22, 23) seem to be non-extensive. However, there is an explicit volume factor, 
N T Nf, in Eq. (24). The resolution is that away from a critical point the variances and 
covariances of the plaquettes scale as 1/V, which is a consequence of the central limit 
theorem. 

Certainly if each plaquette variable could be considered to be fluctuating randomly 
around its mean value then the application of the central limit theorem would be clear. Be- 
fore proceeding, we emphasize that both the plaquette variables defined here are summed 
over all spatial orientations, and hence are invariant under spatial rotations. In the notation 
of Ref. [12], they are projected on the channel. Thus, their covariances are integrals 
over the A± + plaquette correlation function. If plaquette correlations had a finite range, 
then again these terms would be linear in volume if N s were sufficiently large. However, 
if the correlation length associated with plaquettes becomes infinite, then, in the ther- 
modynamic limit, this term would grow faster than the remainder. Consistently, at a second 
order phase transition, where this is expected, C v , as defined in Eq. (18) would scale non- 
trivially with volume according to the critical exponents of the theory. Such behaviour has 
been found in the SO(3) gauge theory [34]. 



2.3 Final expressions 



Expressions for the energy density and the pressure in the usual form are obtained from 
Eq. (13) by multiplying by appropriate powers of N T . In the isotropic (£ = 1) limit and 
for 3+1 dimensions we get 



jr, = 6iV c iV T 4 



Ds-Dr 

9 2 

D x - D T 



(c' s D s + c' T D T ) 
- (c' s D s + c' T D T ) 



+ 6N C N, 



27ra? 



D S + D T 



and 
(25) 



On comparing these expressions with those obtained using the s-favoured scheme [22], 
one can easily see that the new expression for pressure is exactly 1/3 of the old expression 
of the energy density. Since the energy density in the s-favoured scheme comes out to be 
non-negative at all temperatures and on all temporal sizes N T , our new expression for the 
pressure is therefore expected to give non-negative pressure always. The expression for the 
interaction measure 



A 



rp4 



6n c n: 



B(a s ) 
27ra? 



D S + D T 



(26) 
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is same, and also positive, for both the cases. Since both the pressure and the interaction 
measure are non-negative in the t-favoured operator formalism, the energy density must 
also be non-negative. 

Note that A contains B (a s ) as a factor, but this explicit breaking of conformal symmetry 
may be compensated by the vanishing of the factor D s + D T . To determine the coupling 
g 2 , throughout this work, we use the method suggested in Ref. [35], where the one-loop 
order renormalized couplings have been evaluated by using ^-scheme [36] and taking care 
of the scaling violations due to finite lattice spacing errors using the method in Ref. [37]. 

The expressions for £ and a derivatives of F(£, a) in Eq. (22) can be combined by using 
the form of the lattice derivatives in Eq. (8) to get the temperature derivative of F(£, a). 
Finally inserting the derivatives of the coupling (see Eq. 12 and Eq. 21), taking the £ — > 1 
limit, and specializing to d = 3 we get — 



dF 
dT 



B(a s ) 
2ira 2 



[D T - D s ] + 6N c N r N : 



B(a s ) 



2er, 



6N r N T N? 



B(a s ) 
2ira 2 



2-na 2 

+ c' s a ss + c' t <t tt + (c' s + c' T )a s 



(27) 



Proceeding in the similar way as before, in the £ 



dG 
df 



Ds+Dt 

9 2 

- 6N C N T N* 
+ 6N c N T Ng 
+ 6N C N T N* 



c' s D s + 3cJ T D T + c'^D s + 4'£>. 

2(c' T a T , n 



1 limit in d = 3, one obtains — 

B(a s ) 



2ira 2 



[D T - D s 



+ <7 T 



2n H 



+ 



/2 



+ c' T 2 a T ^ T + 2c' s c' T a s 



B(a s ) 
2na 2 



dF 
dT 

+ « 



c' T )a s 



(28) 



For g 



0, i.e. , in the weak-coupling limit, the dominant contribution to all the pla- 
quettes is of order g 2 [38]. Hence, in this limit, Di oc g 2 , and A/T 4 oc g 2 . In the weak- 
coupling limit, therefore, A/T d can be neglected in comparison with e/T d . The scaling 
of Di also implies that ct^ oc g 4 , as a result of which F and its temperature derivative 
are negligible in this limit compared to G and its derivative. Consequently, T — > in this 
limit, resulting in C v /T d — > (d + l)e/T d+1 and C' 2 — > l/d. Note that in any conformal 
invariant theory in d + 1 dimensions one has e = dP, i.e. , C = T = 0, and hence, by Eq. 
(18), one has identical results— C 2 = l/d and C v /T d = {d+ l)e/T d+1 . 



2.4 On the method 



While the expressions in Eq. (25) look different from those in Ref. [22], one may argue 
[39] that standard formulae for change of variables (from the set {£, a T } to {£, a s }) can 
be used to show that both the expressions are identical. However, this conclusion follows 
only if one also demands the values of the couplings g 2 and g 2 to remain the same under 
the change of the scale from a s to a T . As we argue below, this is not true when the weak 
coupling expressions [ Eq. (14)] are used for the couplings. 
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Figure 1. A/T 4 as a function of the bare coupling [3 using a non-perturbative 
(squares) and one-loop order perturbative (pentagons) /3-function, B(a s ). The results 
agree for j3 > 6.55. The plaquette values for iV T = 8 and the values of the non-pertur- 
bative /3-function are taken from Ref. [2]. 

As can be seen form Eq. (14) the Karsch coefficients Ci(£)'s are differences between 
the isotropic and anisotropic couplings. Hence they do not depend on the scale a of the 
isotropic lattice, but only on the parameter which quantifies the difference between the 
isotropic and the anisotropic lattice, i.e. , the anisotropy parameter £. Thus a change of 
scale from a s to a T does not change these Karsch coefficients. In Appendix A we prove 
this explicitly. Given that the Karsch coefficients are same for both the t-favoured and 
the s-favoured schemes, from Eq. [14] it follows that the anisotropic coupling constants 
gi{a,0 are different for the two schemes due to the scale dependence of the isotropic 
coupling constant g(a). Therefore the expressions for e and P are different at finite (but 
small) lattice spacing in the two different approaches. Since the s-favoured and t-favoured 
schemes are different due to the scale dependence of the isotropic coupling constant g(a), 
the difference between the expressions in both the schemes goes as In a, compared to the 
l/o 2 cut-off dependence of the lattice Wilson action. Hence, the difference between the 
two methods is tantamount to modifying the operators. Moreover, for the usual choice of 
scale setting by T = l/N T a T , our approach corresponds to the natural choice of scale in 
Eq. [14]. It is expected that the results from both the methods will match for very large 
temporal lattice size N T . However, as is true with the improvement program in general, 
on small lattices the better operators — t-favoured method in this case — should lead to 
results with lesser artifact errors or alternatively positive pressure at even T < T c . 

While the t-favoured method improves the differential method, leading to positive pres- 
sure, it still requires the use of perturbative couplings. On the other hand, the integral 
method evades them but at the cost of the assumption of homogeneity. For small vol- 
umes used in actual simulations, one may feel reassured by its test in form of agreement 
of results with other methods such as the differential method. Note that the expression for 
A/T 4 is identical for both the integral method and the t-favoured scheme. It depends on 
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the /3-function, B(a s ) in Eq. [15]. A non-perturbatively determined /3-function permits 
the integral method to lead to fully non-perturbative EOS. However, one usually fits a phe- 
nomenological ansatz to extract it from a range of couplings 6/g 2 with their associated 
systematic uncertainties. The differential method could also employ such a /3-function but 
for internal consistency we require that the Karsch coefficients and B{a s ) be obtained at 
the same order, i.e. at one-loop order in the present state of the art. 

The two methods must agree if one uses sufficiently small lattice spacings, viz. when 
the use of perturbative couplings is justified in the differential method computation and on 
large enough volumes. A a comparison between the values of A/T 4 extracted for a given 
N T using the two approaches would reveal at what T the two methods become close to 
each other. Using asymptotic scaling, one could also then find the minimum value of N T 
required for the same level of agreement as a function of T. Such a comparison is shown 
in Figure 1, which demonstrates that a bare coupling of (3 > 6.55 should suffice to give an 
agreement between the t-favoured scheme and the integral method. For (3 < 6.55 use of 
one-loop order perturbative Karsch coefficients may give rise to some systematic effects. A 
comparison with the non-perturbatively determined Karsch coefficients [26,40] shows that 
the difference between the perturbative and non-perturbative values are significant. For 
example, while at around (3 = 6.55 the one-loop order perturbative and non-perturbative 
c- differ by ~ 20%, around (3 = 6 this difference increases to <~ 80%. 

In the present work we show that within the framework of differential method it possible 
to get a positive pressure for all temperatures if one uses the better operators of the t- 
favoured scheme. This is so in spite of the use of one-loop order perturbative Karsch 
coefficients. However, the use of one-loop order perturbative Karsch coefficients [33,3] 
may give some systematic effects if the lattice spacing is not small enough. 

3. Simulations and results 

Our simulations have been performed using the Cabbibo-Marinari pseudo-heatbath algo- 
rithm with Kennedy-Pendleton updating of three SU(2) subgroups on each sweep. Pla- 
quettes were measured on each sweep. For each simulation we discarded around 5000 
initial sweeps for thermalization. We found that the maximum value for the integrated au- 
tocorrelation time for the plaquettes is about 12 sweeps for the T = run at (3 = 6 and 
the minimum was 3 sweeps for the T = 3T C run for iV r = 12. Table 1 lists the details of 
these runs. All errors were calculated by the jack-knife method, where the length of each 
deleted block was chosen to be at least six times the maximum integrated autocorrelation 
time of all the simulations used for that calculation. 

In Ref. [41] it was shown that, at sufficiently high temperature, finite size effects are 
under control if one chooses N s = (T/T C )N T +2 for the asymmetric (N T x N^) lattice. We 
have chosen the sizes of the lattices used at finite T based on this investigation. Close to T c 
the most stringent constraint on allowed lattice sizes comes from the screening mass 
determined in Ref. [14]. Among the temperature values we investigated, this screening 
mass is smallest at 1.25T C where it is a little more than 2T. The choice of N s = 2N T + 2 
satisfies this constraint sufficiently. If future work pushes closer to T c , then larger values of 
N s need to be used in view of the further decrease in the screening mass. At T = 
the constraints are simpler because glueball masses are larger, and also smoother functions 
of (3. For the symmetric (Nf) lattices we have chosen N s — 22 as the minimum lattice 
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Table 1. The coupling (f3), lattice sizes (N T x AT 3 ), statistics and symmetric lattice 
sizes (Nf ) are given for each temperature. Statistics means number of sweeps used for 
measurement of plaquettes after discarding for thermalization. 



J- 1 J-c 


P 




Asymmetric Lattice 




Symmetric Lattice 






size 


stat. 


size 


stat. 




6.0000 


8 x 18 3 


1565000 


22 


253000 


u.y 


O.l 3UU 


1U X 22 


/ZD UUU 


22 






6.2650 


12 x 26 3 


c a a r\r\r\ 

504000 


26 


IP/ AAA 

256000 




6.1250 


8 x 18 3 


1164000 


22 


253000 


1.1 


6.2 OU 


11) X 22 


M/UUU 


22 


ZoUUUU 




6.4200 


12 x 26 


1 AAA 

212000 


26 


1 O z' AAA 

136000 




6.2100 


8 x 18^ 


1903000 


22 


O A 1 AAA 

30 l 000 


l.Zj 


O.JoUU 


1U X 22 


o / /UUU 


ZZ 


21 /UUU 




6.5050 


i a . . a / ' 3 

12 x 26 


390000 


26 


"t ,1 AAAA 

240000 




6.3384 


8 x 18 3 


1868000 


22 4 


544000 


1.5 


6.5250 


10 x 22 3 


1 333000 


22 4 


605000 




6.6500 


12 x 26 3 


882000 


26 4 


335000 




6.5500 


8 x 18 3 


2173000 


22 4 


534000 


2.0 


6.7500 


10 x 22 3 


1671000 


22 4 


971000 




6.9000 


12 x 26 3 


1044000 


26 4 


553000 




6.9500 


8 x 26 3 


1300000 


26 4 


433000 


3.0 


7.0500 


10 x 32 3 


563000 


32 4 


148000 




7.2000 


12 x 38 3 


317000 


38 4 


60000 



size and scaled this up with changes in the lattice spacing in accordance with the analysis 
done in Ref. [3]. 

We performed a — > (continuum) extrapolations by linear fits in a 2 oc l/N 2 at all 
temperatures using the three values N T = 8, 10, and 12. In Figure 2(a) we show our data on 
P/T 4 at finite lattice spacings and the continuum extrapolations for different temperatures, 
both above and below T c . We draw attention to the fact that the pressure is positive on each 
of the lattices we have used and also in the a — > limit. It is an interesting piece of lattice 
physics, not relevant to the continuum limit, that the slope of the continuum extrapolation 
changes sign at T c . This is also true of the continuum extrapolation for e/T 4 as shown in 
Figure 2(b). The extrapolation of both P/T 4 and e/T 4 between 1.1T C and 3T C are similar 
to those shown and have therefore been left out of the figure to avoid clutter. 

Similar continuum extrapolations are shown for C v /T 3 and C 2 S in the two panels of 
Figure 3. In all cases, the continuum extrapolations are smooth, and well fitted by a straight 
line in the range of N T used in this study. As mentioned above, it is interesting lattice 
physics to see that for C v /T 3 also, the slope of the continuum extrapolation flips sign at 
T c . This does not happen for C 2 S . Since this is the derivative of the energy density with 
respect to the pressure, the slope of this quantity depends on the slopes of the continuum 
extrapolation of e/T 4 and P/T 4 . 

The results of continuum extrapolations of our measurements are collected in Table 2. 
It is gratifying to note that the pressure and the entropy are not only positive in the full 
temperature range, but also convex functions of T, as required for thermodynamic stability. 

In the various panels of Figure 4 we show a comparison between the continuum ex- 
trapolated results for different quantities obtained using the t-favoured scheme, s-favoured 
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Figure 2. In the panel (a) we show the dependence of P/T 4 on 1/N% for different 
temperature values. In the panel (b) we show the same for e/T 4 . The \-a error band of 
the continuum extrapolations have been indicated by the lines. 




1/N t 2 1/N, 2 



Figure 3. In the panel (a) we show the dependence of CV/T 3 on 1/N% for different 
values of temperature. In the panel (b) we show the same for Cj. The 1-er error band 
of the continuum extrapolations have been indicated by the lines. 



scheme and the integral method. While the results of the t-favoured and the s-favoured 
schemes are obtained from the analysis of our data, the results of the integral method are 
taken form Ref . [2] . 

First we note that unlike the s-favoured differential method, the t-favoured scheme yields 
a positive pressure [Figure 4(a)] at all T. There is apparent agreement between the integral 
and the t-favoured operator method for T > 2T C , both differing from the ideal value by 
about 20%. Only at these temperatures the coupling (3 becomes > 6.55 for all the lattices 
(see Table 1) that has been used to extract the continuum extrapolated values in the t- 
favoured scheme. Hence, from our earlier discussion it is clear that an agreement between 
the two methods is expected to take place at these temperatures. There can be several 
causes for the difference between these two methods closer to T c — (i) The use of one-loop 
order perturbative Karsch coefficients in the t-favoured scheme is probably the primary 
cause for this difference. Use of larger lattices (i.e. larger (3) or inclusion of the effects of 
higher order loops in the Karsch coefficients is expected to improve the agreement, (ii) 
Another possible source of disagreement is that the results for the integral method shown 
here were obtained on coarser lattices [2] than the ones used in this study, (iii) The integral 
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Figure 4. We show comparisons between the continuum extrapolated results of dif- 
ferent thermodynamic quantities for t-favoured scheme (boxes), the s-favoured scheme 
(triangles) and the integral method (line). In panel (d) we show the continuum extrap- 
olated values of the conformal measure C (boxes). In panel (f) we show a comparison 
between our continuum extrapolated results for CV /T 3 (open boxes) and that of 4e/T 4 
(filled boxes). The data for the integral method has been taken form Ref. [2], 



method assumes that the pressure below some /3q, corresponding to some temperature 
T < T c , is zero. By changing (3q one can change the integral method pressure by a 
temperature independent constant. This may restore the agreement close to T c , although in 
that case the agreement at the high-T region may get spoiled, (iv) Also different schemes 
have been used to define the renormalized coupling in the two cases. This can also make 
some contribution to the different results of the two methods. 

Correspondingly, the energy density is harder near T c , showing a significantly lessened 
tendency to bend down. This could indicate a difference in the latent heat determined by 
the two methods. We shall return to this quantity in the future. The entropy density is 
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E/P — B- 
S/T 3 + TE/3T 4 — e- 



(b) 



2 2.5 
T/T„ 



3 3.5 



Figure 5. In panel (a) we show the temperature dependence of the contribution of one 
of the covariance terms in CV /T 3 . In panel (b) we show the the individual contribution 
of the two factors in Eq. (18) for CV /T 3 . See the text for a detailed discussion. 



shown in Figure 4(c). Since this is a derived quantity (see Eq. 5), it has similar features as 
those of P/T 4 ande/T 4 . 

The generation of a scale and the consequent breaking of conformal invariance at short 
distances, of the order of a, in QCD is, of course, quantified by the /3-function of QCD. 
It has been argued in Ref. [3], that the conformal measure, C = A/e, parametrizes the 
departure from the conformal invariance at the distance scale of order 1/T. In Figure 
4(d) we plot C. It is clear that at high temperature, 2-3T c , conformal invariance is better 
respected in the finite temperature effective long-distance theory. Closer to T c conformal 
symmetry is badly broken even in the thermal effective theory. This is consistent with the 
existence of many mass scales in the theory as found in Ref. [14-16]. It is interesting 
to note that the t-favoured scheme yields marginally smaller values of C than the integral 
method. Note also the peak in C just above T c ; this is the reflection of a similar peak in A. 

Figure 4(e) shows the continuum extrapolated results for . At temperatures of 2T C 
and above, the speed of sound is consistent with the ideal gas value within 95% confidence 
limits. It is seen that Cl decreases dramatically near T c . Below T c there is again a rise in 
Cg , the numerical values being very close 10% below and above T c . In future we plan to 
explore in greater detail the region in between. 

The behaviour of CV/T 3 , shown in Figure 4(f), is the most interesting. At 2T C and above 
it disagrees strongly with the ideal gas value, but is quite consistent with the prediction in 
conformal theories that C v /T 3 = 4e/T 4 . Closer to T c , however, even this simplification 
vanishes. The specific heat peaks at T c , consistent with the observation of Refs. [12,42] 
that there is a light mode (the thermal scalar, called the A^ + ) in the vicinity of T c . Below 
T c the specific heat is very small. 

In view of the rise in CV/T 3 near T c , we studied the contributions of the terms contain- 
ing different covariances of the plaquettes. As can be seen from the Eqs. (27, 28), among 
all the terms containing covariances, the term (a ss + a TT — 2a ST )/g A will have the largest 
contribution to C v /T 3 . All the other terms containing the covariances are multiplied ei- 
ther by one of the c' i7 or by B(a s ) I2-ko? s and hence become at least one order of magnitude 
smaller than this term. 

In Figure 5(a) we show the contribution of the above term, as a function of T in the 
continuum limit. It peaks near T c , consistent with the decrease of the + screening 
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Table 2. Continuum values of some quantities at all temperatures we have explored. 
The numbers in brackets are the error on the least significant digit. For the convenience 
of the readers here we also list the numerical values of these quantities for an ideal gas — 
e/T 4 « 5.26, P/T 4 « 1.75, s/T 3 « 7.02, CV/T 3 « 21.06 and C 2 = 1/3. The 
value of the 't Hooft coupling g 2 N c is computed at the scale 2ttT using the Tc/A-^ 
quoted in Ref. [32]. 



T/T c 


9 Z N C 


e/T 4 


P/T 4 


s/T 3 


C v /T 3 


c'l 


0.9 


11.5(3) 


1.09(4) 


0.14(1) 


1.23(5) 


8.0(5) 


0.162(7) 


1.1 


10.4(2) 


4.31(9) 


0.49(1) 


4.80(6) 


26(2) 


0.18(1) 


1.25 


9.8(2) 


4.6(1) 


0.82(2) 


5.4(1) 


25(1) 


0.21(1) 


1.5 


9.0(1) 


4.5(1) 


1.06(4) 


5.6(2) 


22.8(7) 


0.25(1) 


2.0 


8.1(1) 


4.4(1) 


1.26(4) 


5.7(2) 


17.9(7) 


0.31(1) 


3.0 


7.0(1) 


4.4(1) 


1.37(3) 


5.8(1) 


17.9(8) 


0.32(1) 



mass mentioned earlier. Since the lattices that we used are significantly larger than this 
correlation length, we are in the correct regime of volumes where the central limit theorem 
holds for the fluctuations of the plaquettes. The contribution of this term is very small: 
comparable to the errors in C v . The origin of the peak in C v therefore lies elsewhere. 
In Figure 5(b) we separately plot the two factors, e/P and s/T 3 + Te/3T A , in the the 
expression for C v in Eq. 18. The factor s/T 3 +Te/3T 4 is smooth in the whole temperature 
range, and it is the first factor, e/P, which has a peak near T c . Rewriting this as 3/(1 — C), 
we can recognize that the peak in C v is related to that in A. 



4. Discussion 

In this paper we have proposed a modification, viz. the t-favoured scheme, of the differ- 
ential method for the computation of the QCD equation of state. We have shown that this 
improvement gives positive pressure for all temperatures and iV T used, even when the older 
s-favoured differential method [22] gives negative pressure. Note that this is so in spite of 
the use of the same one-loop order perturbative values for the couplings in both cases. Us- 
ing the t-favoured differential method and by extrapolating to the a — > (continuum) limit 
we obtain the energy density and pressure for a pure gluonic theory in the temperature 
range0.9 < T/T c < 3. These differ from their respective ideal gas values by about 20% at 
3T C , and by much more as one approaches T c . On comparing our results with those of the 
integral method [2], we found that ours are larger for T < 2T C . The primary reason behind 
this disagreement seems to be our use of perturbative couplings. Hence the agreement be- 
tween the t-favoured scheme and the integral method is expected to improve by going to 
larger temporal lattice sizes or equivalently to smaller lattice spacings. 

We have also extended the t-favoured scheme to compute the continuum extrapolated 
results of the specific heat at constant volume and the speed of sound. We found that C v 
peaks near T c where, in addition, C s becomes small. Our results are collected together 
in Table 2 and Figure 4. The most robust quantity on the equation of state in all lattice 
computations is A, and the most interesting (and also stable) feature seen to date is the 
peak in A just above T c . Apart from influencing the EOS, it manifests itself as a peak 
in C v . Since C v could be directly measurable through energy or effective temperature 
fluctuations in heavy-ion collisions, understanding A should be one of the prime goals of 
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3T e 2T e 1 .5Tc 1 ,25T C 




Figure 6. In panel (a) we compare the pressures obtained by t-favoured method 
(boxes), integral method (dotted line) and the g 6 ln(l/<?) order perturbative expansion 
(solid line) of The data for the integral method and the perturbative expansion are taken 
from Ref. [2] and Ref. [40] respectively. The values of the T/A w in Ref. [40] has 
been converted to T/T c using the T c /Aju quoted in Ref. [32]. In panel (b) we show 
the deviation of s/so from 3/4 (boxes) as a function of the 't Hooft coupling. We also 
show the prediction of Eq. (1) (solid line). 

theory. Unfortunately, it seems that at present no tool other than lattice computations are 
available for this task. Even models of this important and stable phenomenon are lacking. 

In view of the fact the perturbation theory fails to reproduce the lattice data on EOS, 
specially close to T c , it is interesting to compare our t-favoured scheme results with that of 
the perturbation theory. In Figure 6(a) we compare the pressure obtained in the t-favoured 
method with that from a dimensionally reduced theory, matched with the 4-d theory per- 
turbatively up to order g 6 ln(l/<?) [43]. WritingPss for the ideal gas (Stefan-Boltzmann) 
value of the pressure, the ratio for P/Psb found in the dimensionally reduced theory [43] 
has an undetermined adjustable constant, c. The pressure determined through dimensional 
reduction agrees with our results almost all the way down to T c , for that value of the con- 
stant (c = 0.7) for which it matches with the integral method in the high temperature range. 
In future it would be interesting to check whether an equally good description is available 
in this approach for the full entropy. This would be a non-trivial extension because pertur- 
bation theory misses A completely. The question, therefore, addresses the non-perturbative 
dynamics of the dimensionally reduced theory. 

The strong coupling result in Eq. (1) of Ref. [17] can be compared with our data on the 
entropy density, s/T 3 . This has to be done in an appropriate window of T where the 't 
Hooft coupling g 2 N c is large and C is small. The strong coupling series is an expansion in 
{g 2 N c )~ 1 ' 2 . For N = 4 SYM, the first term vanishes due to a delicate cancellation and the 
series starts with the (g 2 N c )~ 3 ' 2 term [17]. When some of the supersymmetry is broken, 
this cancellation need not occur and the series could start with a term in (g 2 N c )^ 1 ^ 2 . 
Needless to say, the theory we are studying here, pure QCD, lacks supersymmetry. In 
Figure 6(b) we show the deviation of s/so from 3/4 as a function of the 't Hooft coupling 
(s and g 2 N c are listed in Table 2). Also shown is the prediction of Eq. (1). Comparison 
of our data with the latter shows that the AdS/CFT based theory agrees with our data for 
g 2 N c < 9, or in other words for C < 0.3. As a partial summary of our results, we show 
the equation of state in Figure 7 in the form of a plot of P/T 4 against e/T 4 , useful for 
hydrodynamics. In this plot, the ideal gas for fixed number of colours is represented by a 
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Figure 7. The equation of state of QCD matter. The diagonal line denotes possible 
EOS for theories with conformal symmetry. The circle on the diagonal denotes the 
ideal gluon gas, whose EOS in this form is temperature independent. The ellipses de- 
note 66% error bounds on the measured EOS (see Ref. [41]). The ratio of the axes is a 
measure of the covariance in the measurements of e/T 4 and P/T 4 . The wedges pierc- 
ing these ellipses have average slope , and the opening half-angle of these wedges 
indicate the error in . 

single point which is independent of T, and theories with conformal symmetry by the line 
e = 3P. Pure gauge QCD lies close to the conformal line at high temperature, as shown, 
but deviates strongly nearer T c . 

The slope of the wedges piercing the ellipses indicates the speed of sound — when these 
are parallel to the conformal line then C 2 = 1/3. This is clearly the case at high temper- 
ature. However, there is an increasing flattening of the axis, denoting a drop in C 2 as one 
approaches T c . Note that the slope of the curve joining the middle points of the ellipses 
does not give C 2 , since the plot is of e/T 4 against P/T . In a plot of e against P, it would 
have been correct to assume that the slope gives C 2 . 

Two other physically important effects can be read off the figure. First, the softening 
of the equation of state just above T c is shown by the rapid drop in pressure at roughly 
constant e/T 4 . Second, a large latent heat is indicated by the jump between the last two 
points, at almost the same pressure but very different energy densities. 

A final piece of physics can be deduced from the fact that the low temperature phase 
shows a very small P/T 4 at a significantly large value of e/T 4 > 1 just below T c . This is 
an indication that there are very massive modes in the hadron gas which contribute large 
amounts to e without contributing to P. The small value of C v /T 3 at the same T also 
indicates that the energy required to excite the next state is rather large. We have mentioned 
already that the observations just above T c are compatible with the known spectrum of 
excitations in pure gauge QCD [14]. 
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APPENDIX A: Discussion on the Karsch coefficients 

The Karsch coefficients (a) are differences between the anisotropic and isotropic lattice 
couplings and hence do not depend on the scale a of the isotropic lattice, but only on the 
the anisotropic parameter £. One can see this directly from the derivations in Ref. [33], 
where these have been evaluated up to one-loop order in the perturbation theory. For any 
arbitrary £ ^ 1, all integrals contributing in the effective action S e ff(a,^), mentioned 
in Eq. (2.22) of Ref. [33], are independent of the scale the a. The dependence of a are 
only encoded implicitly in the couplings g^ 2 (a,S > ). Hence S e ff(a,t;) of Eq. (2.22) of 
Ref. [33] is equally valid for a = a T . The values of the Karsch coefficients have been 
evaluated by imposing the condition AS e jf — S e ff(a,^) — S e ff(a, 1) = 0, which is 
again independent of the scale a. Hence the one-loop order Karsch coefficients for both 
the case a — a s (s— favoured scheme) and a — a T (t— favoured scheme) are the same. 

Nevertheless, we derive this equality explicitly in the following. Let us assume that the 
one-loop order perturbative expansions for gf's, around the isotropic lattice coupling g, 
have the following forms 



ffi" V.O = 9~\a.) + Ci{0 + 0[g 2 (a s )}, 
9^ 2 (a T ,0=g- 2 (a T ) + a i (0+O[g 2 (a T )}. 



and 



(Al) 



Our claim is that [dci a = [den (£) a ■ In order to prove it we make a Taylor 

series expansion of g i (a S7 £ l ) around a s = a T , at any fixed £ ^ 1 



g t 2 (a s ,£) =g\ 2 (a T ,t) + J2 



(a s - a T ) r 



n=l 



dx r - 



(A2) 



A £ derivative at constant a s , on Eq. (A2) yields 



a n ( 1\ d 



(A3) 



While [dg(a a )/dZ] a> = 0, [dg(a T )/d$ as = [0ff(a a /O/0£] o . ? 0, from Eq. (Al) it 
follows that 



d9i 2 (ar,£) 



9£ 



dg 2 {a T ) 



d_ 

9? 



g 2 (a s ) + J2 



(a T -a s ) n d n g- 2 (a s ) 



n=l 



da* 



+ 
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n=l s vs 



9a" 



9? 



(A4) 



Substituting Eq. (A4) in Eq. (A3) and using relations in Eq. (Al) to calculate the various 
derivatives one obtains 



na: 



1 



71-1 «"„-2 



da 1 } 



na s 

n=l ' 



+ E- 



71=1 



1 - 



-r 


d n g- 2 (a T ) 






da 7 } 






\d n g- 2 (a T ) 











(A5) 



Finally, taking the £ — > 1 limit, ;'.e. setting a s = a r , one has 



9£ 



(A6) 



A variable transformation from {a s ,£} to {a r ,£} gives ^ (d/d£) a 
a T (d/da T )^. Using it on Eq. (A6) one conclusively proves that 



C(3/9?) OT 



9c: 



(A7) 
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